High accuracy simulations of Kerr tails: coordinate dependence and higher multipoles 
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We investigate the late time behavior of a scalar field on a fixed Kerr background using a 2 + 1 
dimensional pseudospectral evolution code. We compare evolutions of pure axisymmetric multipoles 
in both Kerr-Schild and Boyer-Lindquist coordinates. We find that the late-time power-law decay 
rate depends upon the slicing of the background, confirming previous theoretical predictions for 
those decay rates. The accuracy of the numerical evolutions is sufficient to decide unambiguously 
between competing claims in the literature. 

PACS numbers: 04.25.D-,04.30.Nk,04.70.Bw,02.70.Hm 



I. INTRODUCTION 

Following the original work of Price [H , the late time 
behavior of fields propagating on a Schwarzschild black 
hole has been studied for decades. By now there is 
a rather complete theoretical picture of it, with com- 
pelling supporting numerical evidence. Since in the 
Schwarzschild case the background geometry is spheri- 
cally symmetric, perturbations can be expanded in spher- 
ical harmonics. Solutions with different multipole num- 
bers (£, to) evolve independently. For initial data that is 
not time symmetric, each of these perturbations decays 
at late times as ^-(2^+3)^ regardless of the type of pertur- 
bation (scalar, gravitational, electromagnetic, etc). 

In contrast, our understanding of the late-time decay 
on a Kerr background is still incomplete. In fact, it has 
been the subject of controversy over the past decade. 

Since in the spinning case the background is not spheri- 
cally symmetric, additional multipoles are generated dur- 
ing the evolution of an initially pure multipole. The first 
analytical calculations for the decay rate of these dynam- 
ically generated multipoles were done by Barack and Ori 
[1] and Hod [1]. According to their calculations, the in- 
teractions between different multipoles and the black hole 
angular momentum change the decay rate of each of those 
multipoles compared to the Schwarzschild case. Further- 
more, those calculations also predict that in the Kerr case 
the late time decay rate is not universal, in the sense that 
it depends on the type of perturbations. 

For example. Hod Q predicted that a scalar perturba- 
tion with an initial pure multipole structure with indices 
(£, to) on a Kerr background in Boyer-Lindquist coordi- 
nates would be dominated at late times by the following 
multipole component and decay: 

* oc y(A™)t-(2^+3) if ^ === m or ^ = TO+ 1, (1) 

^ ^ if ^ „ ^ > 2 (evCu) , (2) 

* oc y(fc™+i,™)t-(^+™+2) iii^m>2 (odd). (3) 

In contrast, Burko and Khanna Q argued in favor of a 
"simple picture" , in which at late times the decay would 



also be dominated by the lowest multipole that can be 
generated by mode mixing during the evolution, but that 
such a mode decays in the same way that it would on a 
Schwarzschild background. In other words, the authors 
claimed that the details of the multipole interactions do 
not affect the late time decay rate. Therefore, according 
to this simple picture the late time decay for a massless 
scalar field in a spinning Kerr background should simply 
be i-(2^mi„+3)^ where ^min is the lowest multipole that can 
be generated during the evolution. For example, in the 
axisymmetric case an initial perturbation that is sym- 
metric about the equator would produce even € modes 
by mixing, with taun — 0. An antisymmetric initial per- 
turbation would have imXn = 1- 

The simplest case in which these predictions differ is 
the axisymmetric solution of the wave equation with 
£ = 4 initial data. The lowest multipole that can be 
generated during evolution is the monopole term. The 
simple picture argument predicts for it a decay as in 
Schwarzschild, t^^ , while Hod's calculations predict in- 
stead a decay of . Early numerical studies of the 
wave equation on a Kerr background in Boyer-Lindquist 
coordinates by Krivan Q found an approximate decay 
of which was interpreted by some as an approxi- 

mate verification of Hod's prediction. Others suggested 
that Krivan's results were plagued by numerical prob- 
lems and that Krivan's time decay of t~^^^ should not be 
interpreted as supporting Hod's prediction, but as a nu- 
merical artifact, a transient, or both Krivan in fact, 
had reported problems with the angular differentiation in 
the code used in 0]. In particular the decay rate found 
in does not appear to get closer to with increasing 
angular resolution but to t~^'^ , leaving room for ques- 
tioning the interpretation of the result as a confirmation 
of Hod's prediction. 

In order to settle this issue, Burko and Khanna Q per- 
formed simulations of the wave equation on a Kerr back- 
ground using "ingoing Kerr" coordinates and obtained 
a late time decay very close to , concluding that Kri- 
van's results were indeed misleading because of numerical 
artifacts. 
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It turns out, however, that the simulations of Burko 
and Khanna and Krivan do not correspond to the same 
physical scenario. In the spinning case, the characteriza- 
tion of the initial data in terms of a spherical harmonic 
decomposition is not unique, since there is no preferred 
sphere with respect to which such a decomposition should 
be done. In particular, a pure £ = A initial perturbation 
in ingoing Kerr coordinates does not correspond to pure 
£ = 4 initial data in Boyer-Lindquist coordinates, and 
vice-versa. Not only are the spatial coordinates differ- 
ent, but also the slicing of spacetime. In fact, in a more 
recent analytic calculation, Poisson explicitly showed 
that for perturbations of flat spacetime (where, in partic- 
ular, the contribution due to the spin is treated only to 
leading order) the decay does depend on the choice of co- 
ordinates. Furthermore, Poisson was able to obtain both 
Hod's and the simple prediction within this linearized 
scenario, depending on which coordinates are used. 

In Ref. Scheel et al. presented high accuracy evolu- 
tions of the wave equation on a Kerr background using 
pure multipolc initial data in Kcrr-Schild coordinates up 
to £ = 4. Their numerical results confirmed the predic- 
tions of the simple picture with high accuracy; in par- 
ticular, they were consistent with those of Burko and 
Khanna. 

More recently, Gleiser, Price and PuUin Q have pre- 
sented numerical results for linear scalar perturbations 
of a Schwarzschild black hole where the effect of the spin 
is treated perturbatively. This is done through a hierar- 
chy of one-dimensional evolution equations, which take 
into account increasingly higher order corrections due to 
the spin, in what would correspond to Boyer-Lindquist 
coordinates. Their results show that if one starts with 
a pure £ = 4 mode in those coordinates and solves the 
corresponding hierarchy of equations to study the late 
time behavior of the monopole term (which is the one 
that dominates at late times) at the leading order in the 
angular momentum, it decays as in Hod's prediction (i.e., 
as i~^), not as in the simple picture prediction. The ad- 
vantage of this approach is that the resulting equations 
can be solved by a simple numerical code, since they are 
one-dimensional in space. The disadvantage is that the 
treatment is perturbative in the black hole angular mo- 
mentum, leaving open the question of whether the full 
non-linear dependence on the angular momentum would 
change the results. 

The goal of this paper is to study the late time behav- 
ior of fields in a Kerr spacetime by numerically evolving 
a scalar field in axisymmetry, including the full depen- 
dence of the equations on the spin, and unambiguously 
determine what the decay is, and whether it depends on 
the coordinates used. For that purpose we numerically 
evolve a scalar field on a Kerr background, using both 
Kerr-Schild and Boyer-Lindquist coordinates. Our sim- 
ulations are of enough accuracy so as to rule out any 
possibility of numerical artifacts. We discuss the case 
that has been under dispute, which is the simplest one in 
which Hod's and the simple prediction disagree, namely 



pure 1 = 4 initial data. We also go beyond it and analyze 
the ^ = 5, 6, 7, 8 initial data cases. 

Our results for the Kcrr-Schild case extend those of 
Scheel et al. and Burko and Khanna to higher accu- 
racy and higher multipolc initial data. In particular, we 
confirm that the decay in those coordinates is governed 
by the predictions of the "simple picture". For Boyer- 
Lindquist coordinates our results extend those of Krivan 
and of Gleiser, Pullin and Price, again to higher accu- 
racy and higher multipolc initial data. Consistent with 
the conclusions of those references, we unambiguously 
find that the decay rate in these coordinates is the one 
predicted by Hod. 

The organization of this paper is as follows. In Sec. |TT] 
we describe the method we use to numerically evolve the 
scalar field. In Sec. IIIII we present our results for pure £ 
initial data for both Kerr-Schild and Boyer-Lindquist co- 
ordinates. Finally, we present our conclusions in Sec. lIVI 



II. THE NUMERICAL METHOD 

We evolve axisymmetric solutions of the wave equation 
V'^Va^' = on a fixed Kerr black hole background in 
both Kerr-Schild and Boyer-Lindquist coordinates. In 
both cases we cast and numerically solve the equations 
as a first order in space and time system. For example, 
in the Boyer-Lindquist case we write the equations as 



(r^ + a^)^ 
, A(l-y2) 



E2 



5„* 



2rA 



2yA 

S2 ^ 



where M, a are the mass and spin of the black hole, re- 
spectively. 



-2Mr + a'^, 



S = [{r' + ay-a'A{l-y')] 



1/2 



y 



— cos y. 



and where we have introduced the Kerr-tortoise coordi- 
nate r*, defined by 



dr 



Here ^f^, , '^t , and are auxiliary variables introduced 
to cast the system in first order form; at the continuum 
they satisfy ^Pr. = = 9t^', and = dy"^. 

We solve the resulting equations using a special pur- 
pose two-dimensional code written for this project. A re- 
quirement for the numerical solution is that it be demon- 
strably accurate enough that there is no doubt about the 
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results. Since the solution is smooth, a spectral method 
should be optimal in terms of efficiency for high accu- 
racy. Accordingly, we use a pseudo-spectral collocation 
(PSC) method in space and the method of lines to evolve 
in time. A by-product of using a spectral method is that 
we avoid all difficulties in handling the polar singularities 
in spherical coordinates. 

In the simulations shown below the domain in the ra- 
dial direction is partitioned into blocks, each of length 
lOM (in either Boyer-Lindquist-tortoise or Kerr-Schild 
coordinates). On each block the dependence of the so- 
lution in the angular and radial directions is expanded 
in spherical harmonic and Chebyshev polynomials, re- 
spectively, using Gauss-Lobatto collocation points. The 
solution is advanced in time at each of these points us- 
ing a fourth-order Rungc-Kutta method. Information 
at the interfaces of the different blocks is communicated 
through characteristic variables using a penalty method, 
as described in [§, [l^] . 

We use initial data of the form 



1=4 Kerr-Schild initial data and evolution 
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Figure 1: Decay of an initial £ = 4 mode in Kerr-Schild coor- 
dinates. 



III. RESULTS 



«'(t = 0) = 0, 

where Tq = 20M and a = AM. In these expressions as 
well as in the results shown in the next section, r and 
t are the Boyer-Lindquist or Kerr-Schild radial and time 
coordinates, respectively, depending on the equations be- 
ing solved. 

We set the angular momentum to a = 0.5A/. In the 
Kerr-Schild case the black hole singularity is excised by 
placing a purely outflow inner boundary at 1.8A/, while 
in the Boyer-Lindquist case the inner boundary is placed 
at = — 40M (which in Boyer-Lindquist radius corre- 
sponds to a distance of ~ 10^^ M from the event horizon) 
and we set the incoming modes to zero. In both cases 
the outer boundary is placed far away enough so that 
the results here shown are causally disconnected from 
the type of boundary conditions there imposed (incom- 
ing modes set to zero). In more detail, if wc evolve for 
a total time T, we typically place the outer boundary at 
r = T/2 + lOOAf -|- rjB, where rjB is the radial location 
of the inner boundary. 

The number of collocation points in the radial and an- 
gular directions per domain is denoted by nr,ni, respec- 
tively and; unless otherwise stated, the time step is kept 
fixed at At — 0.025A/. Obtaining fast convergence as 
we increase the number of collocation points shows that 
the time step is small enough that the errors due to the 
time integration are smaller than those associated with 
the spatial dimensions. 

In order to be able to follow the solution for long 
enough periods of time, especially for the higher mul- 
tipole initial data cases (which decay faster) we use 
quadruple precision. 



Hod's predictions and the simple one coincide for £ = 
0,1,2,3 initial data. We have checked that our simu- 
lations reproduce the expected decay for each of those 
values of £, using both Boyer-Lindquist and Kerr-Schild 
coordinates. In particular, we have found that £ — and 
£ = 1 modes already present in the initial data have a late 
time decay of and t~^, respectively. As we will show 
below, and as predicted by Hod and Poisson, the decay 
of these multipoles is different when they are dynamically 
generated in a Boyer-Lindquist background. 

In the following we report our results for higher mul- 
tipolc initial data, for which Hod's prediction and the 
simple one disagree. Wc show our results for observers 
at r = 21.8M in the Kerr-Schild case and r» = 20 Af in 
the Boyer-Lindquist one, though similar results hold for 
other observer locations. 



A. The £ = 4 case 

Figure [1] shows the decay of the scalar field versus time 
for a fixed observer radius and £ = 4 initial data and 
evolution in Kerr-Schild coordinates. Shown arc the first 
few {£ = 0,2,4) even multipole components of the so- 
lution [the odd components stay at quadruple precision 
roundoff values (~ 10"'^^) at all times]. Initially only the 
£ = 4 mode is present, but additional modes are gener- 
ated during evolution; at late times the monopole term 
dominates. 

Figure [2] shows the local power index, defined as 

LPI{£,t) = -t*fV^^^^ (where is the projection 
of the solution to its i?-multipole component), for the 
monopole term shown in FiglT] If the monopole term de- 
cays as (X t^^ at late times, then one should obtain 
that LPI — > /X. The LPI for the monopole is approaching 
—3 at late times in our simulations, as predicted by the 
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1=4 Kerr-Schild initial data and evolution 



1=4 Boyer-Lindquist initial data and evolution 
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Figure 2: Local power index for the monopole term of the 
previous figure. 
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Figure 4: Local power index for the monopole term of the 
previous figure. 



1=4 Boyer-Lindquist initial data and evolution 
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Figure 3: Decay of an initial I = A mode in Boyer-Lindquist 
coordinates. 



simple picture. 

Figure [3] shows the decay for the evolution of £ = 4 
initial data in Boyer-Lindquist coordinates. The slowest 
decaying multipolc is again the monopole term, but in 
this case we find that it decays as t^^, as predicted by 
Hod and shown in Fig. |4l 

The results shown above were obtained with = 45 
and rii = 4. To give an idea of the errors in these 
simulations, Fig. [5] shows the differences in the LPI for 
the monopole term between different spatial resolutions 
{rir = 30,35,38,40) and the highest one (n^ = 45) as 
a function of time, keeping ~ A fixed. Similarly, the 
figure also shows the differences in the LPf when keeping 
rir = 45 fixed and increasing n^. At late times the errors 
introduced by using Legendre polynomials, where i 
is the multipole index of the initial data (as opposed to 
using higher order polynomials) is comparable to the er- 
rors in the radial direction. Therefore, unless otherwise 
specified, in the simulations below we use ni Legendre 
polynomials, where the value of £ is the same as in the 



B. The 



5 case 



Figure [S] shows our results for the dipole LPI in evo- 
lutions of ^' = 5 initial data in Kerr-Schild and Boyer- 
Lindquist coordinates. In both cases the dipole term 
eventually dominates, with a decay of t^^ in the Kerr- 
Schild case, as predicted by the "simple picture" , and a 
decay of in Boyer-Lindquist coordinates, as predicted 
by Hod. 



C. The 1 = 6 case 

Figure [7] shows our results for the monopole LPI in 
evolutions of ^ = 6 initial data in Kerr-Schild and Boyer- 
Lindquist coordinates. According to the simple picture 
interpretation, the monopole term should dominate at 
late times, with a decay of t^^ , while Hod's prediction 
for this case is a decay of 



D. The 1 = 7 case 

Figure [8] shows our results for the dipole LPI in evo- 
lutions of £ = 7 initial data in Kerr-Schild and Boyer- 
Lindquist coordinates. In the latter case the solution de- 
cays much faster and we need to use higher resolution in 
the radial direction, with an associated smaller timestep 
for CFL stability. Figure [5] shows our results for both 
71,. = 45, M = 0.025 and rir = 50, At = 0.0125. Accord- 
ing to the simple picture interpretation, the dipole term 
should dominate at late times, with a decay of , while 
Hod's prediction for this case is a decay of . 
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1=4 BL initial data and evolution 



1=5 Kerr-Schild initial data and evolution 
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Figure 5: Errors in the LPI for the monopole term of the 
previous figure. 
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1=5 Boyer-Lindquist initial data and evolution 
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Figure 6: Local power index for the dipole term of ^ = 5 
Kerr-Schild and Boyer-Lindquist initial data and evolution. 



E. The ^ = 8 case 

In the Boyer-Lindquist case the monopole term decays 
even faster for £ = 8 initial data and we start reach- 
ing quadruple precision roundoff values by the time the 
tail regime begins, as seen in Figure [5] As a conse- 
quence, there is a rather short time interval in which 
we can measure the local power index, and the latter is 
not as clean as for the lower multipoles. Figure [TO] shows 
our results for evolutions of ^ = 8 initial data in both 
Kerr-Schild and Boyer-Lindquist coordinates. According 
to the simple picture interpretation, the monopole term 
should dominate at late times, with a decay of t~^, while 
Hod's prediction for this case is a decay of t~^. 



IV. DISCUSSION 

In this paper we have evolved an axisymmetric scalar 
field on a Kerr background using both Boyer-Lindquist 
and Kerr-Schild coordinates, with pure multipolc initial 
data with indices £ = 0, 1, 2 ... 6, 7, 8, and established the 
decay rate at late times for each initial data and choice of 



background coordinates. The numerical values that we 
obtained for those rates confirm the "simple picture" pre- 
diction Q in the Kerr-Schild case, but Hod's prediction 
Q in the Boyer-Lindquist one. Thus, the most impor- 
tant conclusion of this paper is that the observed late- 
time decay of the field depends on the time-slicing of the 
background spacetime. 

The differences between the decay rates that we found 
in our simulations and the above asymptotic {t — > oo) 
predictions are small enough so as to rule out the possi- 
bility of numerical artifacts. Those differences are typi- 
cally less than one percent, except for the £ = 8 initial 
data Boyer-Lindquist case. The dominant factor in these 
differences in all cases appears to be the fact that we 
are evolving for a long but finite time; the numerical er- 
rors in our simulations are several order of magnitudes 
smaller. In the particular case of Boyer-Lindquist evolu- 
tions of £ = 8 initial data, the solution reaches quadruple 
precision roundoff errors while entering the tail regime. 
For that reason we cannot determine the decay rate with 
uncertainties as small as in the £ = ... 7 cases. The 
late time decay rate that we find for the £ — 8 initial 
data Boyer-Lindquist case is roughly between —9.2 and 
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1=6 Kerr-Schild initial data and evolution 



1=7 Kerr-Schild initial data and evolution 
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1=6 Boyer-Lindquist initial data and evolution 
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Figure 7: Local power index for the monopole term of £ — 6 
Kerr-Schild and Boyer-Lindquist initial data and evolution. 



Figure 8: Local power index for the dipole term of ^ = 7 
Kerr-Schild and Boyer-Lindquist initial data and evolution. 



—9.0, to be compared with Hod's prediction 
the simple picture one (—3). 

Our results are in apparent, but not real, contradic- 
tion with those of Scheel et al. Q and recent work by 
Burko and Khanna . In order to clarify the source of 
these apparent contradictions we need to describe parts 
of those references in some detail. 

The 3D code used by Scheel et al. in Ref. 0] only 
allowed for evolutions on a black hole background where 
the singularity is excised from the computational domain 
by placing a purely outflow inner boundary. Since this 
excludes Boyer-Lindquist coordinates, Kerr-Schild ones 
{t, X, y, z) were used in Ref. 0]. In this coordinate system, 
the metric is given by 



9) and coordinate but rather a spheroidal coordinate: 



where 



H = 



4 I 9 9' 



, rx + ay ry — ax z 



and the radial coordinate r is defined not as a spherical 



x^ + ^ 
+ a^ 

Scheel et al. evolved scalar waves on a Kerr-Schild back- 
ground using two sets of initial data, corresponding to 
pure multipoles in spherical and spheroidal coordinates 
on constant-time Kerr-Schild slices. In both cases, they 
found the tail decayed as in the simple picture. Our 
Kerr-Schild evolutions confirm this. However, this does 
not necessarily correspond to the evolution of pure multi- 
pole initial data on a constant-time Boyer-Lindquist slice, 
as seen by a Boyer-Lindquist observer. Hod's prediction 
is for the late time decay rate seen by an observer with 
constant Boyer-Lindquist radius as a function of Boyer- 
Lindquist time, for an evolution of pure multipole ini- 
tial data in Boyer-Lindquist coordinates. This is exactly 
what we have explicitly numerically solved for in this pa- 
per. 

In more recent work (ll| . Burko and Khanna evolved 
pure multipole f = 4 initial data, using Boyer-Lindquist 
coordinates to define the multipole decomposition as well 
as to perform the evolution, obtaining a decay of t~^, 



7 



1=8 Boyer-Lindquist initial data and evolution 
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Figure 9: Decay of the monopole term for £ — 8 Boyer- 
Lindquist initial data and evolution. 



which appears to be consistent with the simple picture 
and in contradiction with our resuhs. The reason for 
this apparent contradiction seems to be caused by the 
type of initial data used in Ref . ^lli] . Rather than giving 
pure multipole initial data to the scalar field and its time 
derivative as done in Hod's original work and in this pa- 
per, in Ref. [ll| pure multipole initial data was given to 
the scalar field and its "momentum" , 



0) 
0) 



0, 



where b = b{r, 9) . The angular dependence of b effec- 
tively corresponds to adding an £ = component to the 
initial data for dt"^. As explained at the beginning of 
Sec. mil the decay rate of an initial £ = mode for the 
pair (^f, i9tVE') according to both Hod's prediction and the 
simple picture is (our simulations, as well as previ- 
ous ones, have in particular confirmed this). What we 
have found in this paper, though, which is also reported 
by Gleiser, Price and Pullin 0], is that the decay rate of 
a dynamically generated monopole in a Boyer-Lindquist 
coordinates is faster. If as initial data for (5", 9(^1') one 
superposes £ = A and £ ~ modes, as effectively done in 
Ref. [ll| , the late time decay rate of each component can 
be considered independently, since the evolution equa- 
tions considered are linear. The evolution of the initial 
£ = 4 component will dynamically generate a monopole 
term that decays as t~^, while the evolution of the initial 
£ = component will be dominated by a decay rate of 
t^^. The results of Ref. [ll| can be understood by the 
fact that the latter will dominate over the former at late 



times. 
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